Imports System.IO
Imports MatrixDataStructures

#Region "PAR file parameters"
Class PARfileInfo
    'Public dataset_name As String
    'Public patient_name As String
    'Public examination_name As String
    'Public protocol_name As String
    'Public exam_time_date As String
    'Public acquisition_nr As Integer
    'Public reconstruction_nr As Integer
    'Public scan_duration_secs As Integer
    'Public max_no_of_cardic_phases As Integer
    'Public max_no_of_echoes As Integer
    Public max_no_of_slices As Integer
    'Public max_no_of_dynamics As Integer
    'Public max_no_of_mixes As Integer
    Public image_pixel_size As Integer
    'Public technique As String
    'Public scanmode As String
    'Public scan_resolution_x_y(2) As Integer
    'Public scan_percentage As Integer
    Public recon_resolution_x_y(2) As Integer
    'Public no_of_averages As Integer
    'Public repetition_time As Single
    Public fov(3) As Single
    Public slice_thickness As Single
    Public slice_gap As Single
    'Public water_fat_shift As Single
    Public angulation_midslice(3) As Single
    Public off_center_midslice(3) As Single
    'Public flow_compensation As Integer
    'Public pre_saturation As Integer
    'Public cardiac_frequency As Integer
    'Public min_rr_interval As Integer
    'Public max_rr_interval As Integer
    'Public phase_encoding_velocity(3) As Single
    'Public mtc As Integer
    'Public spir As Integer
    'Public epi_factor As Integer
    'Public turbo_factor As Integer
    'Public dynamic_scan As Integer
    'Public diffusion As Integer
    'Public diffusion_echo_time As Single
    'Public inversion_delay As Single
End Class
Public Class RECimagesInfo
    'other rec file value
    'these values are assumed to be same for all images in the rec file
    'Unneeded Values are Hashed
    'Public image_flip_angle As Single
    'Public diffusion_b_factor As Single
    'Public trigger_time As Single
    'Public dyn_scan_begin_time As Single
    'Public echo_time As Single
    Public slice_number As Integer
    'Public echo_number As Integer
    'Public dynamic_scan_number As Integer
    'Public cardiac_phase_number As Integer
    Public image_type_mr As Integer
    'Public scanning_sequence As Integer
    Public index_in_rec_file As Integer
    Public rescale_intercept As Single
    Public rescale_slope As Single
    Public scale_slope As Single
    'Public window_center As Integer
    'Public window_width As Integer
    'Public image_angulation(3) As Single
    'Public offcenter(3) As Single
    'Public image_orientation As Integer
    Public slice_orientation As Integer
    'Public fmri_status_indication As Integer
    'Public image_type_ed_ed As Integer
    'Public pixel_spacing(2) As Single
End Class
#End Region

#Region "REC file structures"
Class rec_matrix
    Public pixel_size As Integer
    Friend x As Single
    Friend y As Single
    Friend z As Single
    Public volume(,,) As Single
    Public Sub initialize_structure(ByVal a As Integer, ByVal b As Integer, ByVal c As Integer)
        x = a
        y = b
        z = c
        volume(CInt(x - 1), CInt(y - 1), CInt(z - 1)) = New Single
        Dim i, j, k As Integer
        For i = 0 To CInt(x - 1)
            For j = 0 To CInt(y - 1)
                For k = 0 To CInt(z - 1)
                    volume(CInt(x), CInt(y), CInt(z)) = 0.0F
                Next
            Next
        Next
    End Sub
End Class
#End Region

#Region "PARREC Loading Functions"
Module LoadingFunctions
    Function return_nos_from_string(ByVal s As String, ByVal n As Integer) As Single()
        s = s & "  "
        Dim answer(n - 1) As Single
        Dim string_len As Integer
        string_len = Len(s)
        Dim char_array(string_len) As Char
        char_array = s.ToCharArray
        Dim i As Integer = 0
        Dim last_blank_pos As Integer = 1
        Dim counter As Integer = 0
        While (i < string_len - 1)
            If char_array(i) = " " And char_array(i + 1) = " " Then
                answer(counter) = Single.Parse(Mid(s, last_blank_pos, i - last_blank_pos + 1))
                counter = counter + 1
                last_blank_pos = i + 2
            End If
            i = i + 1
        End While
        Return answer
    End Function
    Public Function single_volume_to_3d(ByVal a As Matrix4DInt16, ByVal slice As Integer) As Matrix3DInt16
        Dim ans As New Matrix3DInt16(a.i, a.j, a.k)

        Dim i, j, k As Integer
        For i = 0 To ans.x - 1
            For j = 0 To ans.y - 1
                For k = 0 To ans.z - 1
                    ans.data(i, j, k) = a.data(i, j, k, slice - 1)
                Next
            Next
        Next
        Return ans
    End Function
    Public Function reshape(ByVal matrix(,,) As Integer, ByVal matrix_x As Long, ByVal matrix_y As Long, ByVal matrix_z As Long, ByVal x As Integer, ByVal y As Integer, ByVal z As Integer, ByVal k As Integer) As Matrix4DInt16
        'RESHAPE(X,M,N) returns the M-by-N matrix whose elements
        'are taken columnwise from X.  An error results if X does
        'not have M*N elements.

        'it works like this:
        'it reads the first column of the first slice
        'then to the 2nd and then 3rd column of the first slice
        'then proceeds to the next slice
        'then to the next volume
        Dim data As New Matrix4DInt16(x, y, z, k)
        Dim temp As Long = matrix_x * matrix_y * matrix_z
        Dim array_list(CInt(temp)) As Integer
        If temp <> x * y * z * k Then
            MsgBox("Invalid call of function")
            Return Nothing
        Else
            'this part reads column wise in a array_list
            Dim a, b, c As Long
            Dim d As Long
            d = 0
            For c = 0 To matrix_z - 1
                For b = 0 To matrix_y - 1
                    For a = 0 To matrix_x - 1
                        If matrix(CInt((a)), CInt((b)), CInt(c)) = -32711 Then MsgBox(d & "asds")
                        array_list(CInt(d)) = matrix(CInt(a), CInt(b), CInt(c))
                        d += 1
                    Next
                Next
            Next

            'this part writes column wise from the array_list
            Dim counter As Long = 0
            For d = 0 To k - 1
                For c = 0 To z - 1
                    For b = 0 To y - 1
                        For a = 0 To x - 1
                            If array_list(CInt(counter)) > 32700 Then
                                array_list(CInt(counter)) = 32700
                                MsgBox("overflow corrected!")
                            End If

                            data.data(CInt(a), CInt(b), CInt(c), CInt(d)) = CShort(array_list(CInt((counter))))
                            If data.data(CInt(a), CInt(b), CInt(c), CInt(d)) = -32711 Then MsgBox(counter & "saaaaaa")
                            counter = counter + 1
                        Next
                    Next
                Next
            Next
            Return data
        End If
    End Function
End Module
#End Region

Public Class SingleOptimizedParRecLoader
    Dim LoadedRawData As SingleVolume
    Public Sub New(ByVal FileName As String)
        LoadedRawData = New SingleVolume
        Dim Dimensions As New Coordinate(0, 0, 0)
        'read the header list
        Dim dti_slice_info(0) As RECimagesInfo
        Dim pixel_size, orient As Integer
        Dim rs, ri As Single 'needed to produce display 
        OpenOptimizedPAR(FileName, dti_slice_info, LoadedRawData.Center, LoadedRawData.OffCenterAngle, Dimensions, LoadedRawData.FOV, pixel_size, rs, ri, orient)

        'origin, dimensions and fov are already set
        'nx ny nz as dimensions

        'load the data into a specific volume
        LoadedRawData.Matrix = OpenOptimizedREC(FileName, Dimensions, pixel_size, rs, ri, orient)
        'return the object
    End Sub
    Public Function GetData() As SingleVolume
        If LoadedRawData Is Nothing Then MsgBox("!")
        Return LoadedRawData
    End Function

    Public Function GetRAWData() As Short(,,)
        Return LoadedRawData.Matrix.data
    End Function
    Sub OpenOptimizedPAR(ByVal fname As String, ByRef recc_info() As RECimagesInfo, ByRef origin As CoordinateD, ByRef originAngle As CoordinateD, ByRef dimensions As Coordinate, ByRef fov As CoordinateD, ByRef pixel_size As Integer, ByRef rs As Single, ByRef ri As Single, ByRef orientation As Integer)
        Dim fs As FileStream 'Declare new stream

        fs = New IO.FileStream(fname.ToString, FileMode.Open)
        Dim r As StreamReader = New StreamReader(fs) 'create reader
        'read data into structure
        Dim par_struct As New PARfileInfo
        Dim temp As String
        Dim temp2() As Single
        r.ReadLine()                            '# === DATA DESCRIPTION FILE ======================================================
        r.ReadLine()                            '#
        r.ReadLine()                            '# CAUTION - Investigational device.
        r.ReadLine()                            '# Limited by Federal Law to investigational use.
        r.ReadLine()                            '#
        temp = r.ReadLine()                     '# Dataset name: e:\intera\TH_E004_3_1
        'par_struct.dataset_name = Mid(temp, 17, Len(temp) - 17 + 1)
        r.ReadLine()                            '#
        temp = r.ReadLine()                            '# CLINICAL TRYOUT             Research image export tool     V3
        If temp.Substring(61, 2) = "V3" Then
            'MsgBox("V3")
            r.ReadLine() '#
            r.ReadLine()                            '# === GENERAL INFORMATION ========================================================
            r.ReadLine()                            '#
            temp = r.ReadLine                       '.    Patient(Name)                      :   TH_E004_11292004()
            'par_struct.patient_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination(Name)                  :   Eureka()
            'par_struct.examination_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Protocol name                      :   T1W/3D/TFE
            'par_struct.protocol_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination date/time              :   2004.11.29 / 15:27:35
            'par_struct.exam_time_date = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Acquisition nr                     :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Reconstruction nr                  :   1
            'par_struct.reconstruction_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan Duration [sec]                :   113
            'par_struct.scan_duration_secs = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Max. number of cardiac phases      :   1
            'par_struct.max_no_of_cardic_phases = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of echoes              :   1
            'par_struct.max_no_of_echoes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of slices/locations    :   60
            par_struct.max_no_of_slices = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of dynamics            :   1
            'par_struct.max_no_of_dynamics = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of mixes               :   1
            'par_struct.max_no_of_mixes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Image pixel size [8 or 16 bits]    :   16
            par_struct.image_pixel_size = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            pixel_size = CInt(par_struct.image_pixel_size)
            temp = r.ReadLine                       '.    Technique                          :   T1TFE()
            'par_struct.technique = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan mode                          :   3D
            'par_struct.scanmode = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan resolution  (x, y)            :   240  150
            'temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 2)
            'par_struct.scan_resolution_x_y(0) = CInt(temp2(0))
            'par_struct.scan_resolution_x_y(1) = CInt(temp2(1))
            temp = r.ReadLine                       '.    Scan percentage                    :   100
            'par_struct.scan_percentage = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Recon resolution (x, y)            :   256  256
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 2)
            par_struct.recon_resolution_x_y(0) = CInt(temp2(0))
            par_struct.recon_resolution_x_y(1) = CInt(temp2(1))
            temp = r.ReadLine                       '.    Number of averages                 :   1
            'par_struct.no_of_averages = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Repetition time [msec]             :   9.31  
            'par_struct.repetition_time = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    FOV (ap,fh,rl) [mm]                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.fov = temp2
            temp = r.ReadLine                       '.    Slice thickness [mm]               :   2.00
            par_struct.slice_thickness = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Slice gap [mm]                     :   0.00
            par_struct.slice_gap = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Water Fat shift [pixels]           :   3.05
            'par_struct.water_fat_shift = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Angulation midslice(ap,fh,rl)[degr]:   0.00  0.00  0.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.angulation_midslice = temp2
            temp = r.ReadLine                       '.    Off Centre midslice(ap,fh,rl) [mm] :   -0.16  45.35  0.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.off_center_midslice = temp2
            temp = r.ReadLine                       '.    Flow compensation <0=no 1=yes> ?   :   0
            'par_struct.flow_compensation = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Presaturation     <0=no 1=yes> ?   :   1
            'par_struct.pre_saturation = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Cardiac frequency                  :   0
            'par_struct.cardiac_frequency = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Min. RR interval                   :   0
            'par_struct.min_rr_interval = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. RR interval                   :   0
            'par_struct.max_rr_interval = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Phase encoding velocity [cm/sec]   :   0.00  0.00  0.00
            'temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            'par_struct.phase_encoding_velocity = temp2
            temp = r.ReadLine                       '.    MTC               <0=no 1=yes> ?   :   0
            'par_struct.mtc = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    SPIR              <0=no 1=yes> ?   :   0
            'par_struct.spir = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    EPI factor        <0,1=no EPI>     :   1
            'par_struct.epi_factor = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    TURBO factor      <0=no turbo>     :   150
            'par_struct.turbo_factor = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Dynamic scan      <0=no 1=yes> ?   :   0
            'par_struct.dynamic_scan = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Diffusion         <0=no 1=yes> ?   :   0
            'par_struct.diffusion = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Diffusion echo time [msec]         :   0.00
            'par_struct.diffusion_echo_time = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Inversion delay [msec]             :   0.00
            'par_struct.inversion_delay = CSng(Mid(temp, 45, Len(temp) - 45 + 1))

            'now load info about the rec images

            r.ReadLine() '#
            r.ReadLine() '# === PIXEL VALUES =============================================================
            r.ReadLine() '#  PV = pixel value in REC file, FP = floating point value, DV = displayed value on console
            r.ReadLine() '#  RS = rescale slope,           RI = rescale intercept,    SS = scale slope
            r.ReadLine() '#  DV = PV * RS + RI             FP = DV / (RS * SS)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION DEFINITION =============================================
            r.ReadLine() '#  The rest of this file contains ONE line per image, this line contains the following information:
            r.ReadLine() '#
            r.ReadLine() '#  slice number      	                       (integer)
            r.ReadLine() '#  echo number                              (integer)
            r.ReadLine() '#  dynamic scan number                      (integer)
            r.ReadLine() '#  cardiac phase number                     (integer)
            r.ReadLine() '#  image_type_mr                            (integer)
            r.ReadLine() '#  scanning sequence                        (integer)
            r.ReadLine() '#  index in REC file (in images)            (integer)
            r.ReadLine() '#  rescale intercept                        (float)
            r.ReadLine() '#  rescale slope                            (float)
            r.ReadLine() '#  scale slope                              (float)
            r.ReadLine() '#  window center                            (integer)
            r.ReadLine() '#  window width                             (integer)
            r.ReadLine() '#  image angulation (ap,fh,rl in degrees )  (3*float)
            r.ReadLine() '#'  image offcentre (ap,fh,rl in mm )        (3*float)
            r.ReadLine() 'r.ReadLine()#'  image_display_orientation                (integer)
            r.ReadLine() '#  slice orientation ( TRA/SAG/COR )        (integer)
            r.ReadLine() '#  fmri_status_indication                   (integer)
            r.ReadLine() '#  image_type_ed_es  (end diast/end syst)   (integer)
            r.ReadLine() '#  pixel spacing (x,y) (in mm)              (2*float)
            r.ReadLine() '#  echo_time                                (float)
            r.ReadLine() '#  dyn_scan_begin_time                      (float)
            r.ReadLine() '#  trigger_time                             (float)
            r.ReadLine() '#  diffusion_b_factor                       (float)
            r.ReadLine() '#  image_flip_angle (in degrees)            (float)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION ==========================================================
            r.ReadLine() '#sl ec dyn ph ty  idx (re)scale             window       angulation        offcentre         info     spacing   echo  dtime ttime diff  flip'
            r.ReadLine() '..

            'fill up the slices combo box and data in the array
            Dim i As Integer = 0
            Dim array_list As New ArrayList
            temp = r.ReadLine.ToString

            Dim start As Integer = 0
            Dim temp3 As String

            Do While temp.Length > 100
                'recc_info(i) = New RECimagesInfo
                'declare array of objects and fill it
                Dim list_data As New RECimagesInfo
                list_data = New RECimagesInfo
                start = 0
                list_data.slice_number = CInt(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)

                list_data.image_type_mr = CInt(ReadNextAfterSpace(temp, start))
                list_data.index_in_rec_file = CInt(ReadNextAfterSpace(temp, start))
                list_data.rescale_intercept = CSng(ReadNextAfterSpace(temp, start))
                list_data.rescale_slope = CSng(ReadNextAfterSpace(temp, start))
                list_data.scale_slope = CSng(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                list_data.slice_orientation = CInt(ReadNextAfterSpace(temp, start))

                array_list.Add(list_data)
                temp = r.ReadLine.ToString
            Loop

            'now that we know the volume, move all the data to the recc_info structure and dispose the array_list
            ReDim recc_info(par_struct.max_no_of_slices)
            Dim a As Integer
            For a = 0 To par_struct.max_no_of_slices * 1 - 1
                recc_info(a) = CType(array_list.Item(a), RECimagesInfo)
            Next
            array_list.Clear()
            fov = New CoordinateD(0, 0, 0)
            fov.x = CShort(CInt(par_struct.fov(2)))
            fov.y = CShort(CInt(par_struct.fov(0)))
            fov.z = CShort(par_struct.fov(1)) '(par_struct.slice_gap + par_struct.slice_thickness) * par_struct.max_no_of_slices
        ElseIf temp.Substring(61, 2) = "V4" Then
            r.ReadLine() '#
            r.ReadLine()                            '# === GENERAL INFORMATION ========================================================
            r.ReadLine()                            '#
            temp = r.ReadLine                       '.    Patient(Name)                      :   TH_E004_11292004()
            'par_struct.patient_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination(Name)                  :   Eureka()
            'par_struct.examination_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Protocol name                      :   T1W/3D/TFE
            'par_struct.protocol_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination date/time              :   2004.11.29 / 15:27:35
            'par_struct.exam_time_date = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Series Type                      :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Acquisition nr                     :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Reconstruction nr                  :   1
            'par_struct.reconstruction_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan Duration [sec]                :   113
            'par_struct.scan_duration_secs = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Max. number of cardiac phases      :   1
            'par_struct.max_no_of_cardic_phases = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of echoes              :   1
            'par_struct.max_no_of_echoes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of slices/locations    :   60
            par_struct.max_no_of_slices = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of dynamics            :   1
            'par_struct.max_no_of_dynamics = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of mixes               :   1
            'par_struct.max_no_of_mixes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            r.ReadLine() 'Patient Position
            r.ReadLine() 'Preparation Direction
            r.ReadLine() ' technique
            r.ReadLine() ' scan resolution
            r.ReadLine() 'scan mode
            r.ReadLine() 'repitition time
            temp = r.ReadLine                       '.    FOV (ap,fh,rl) [mm]                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.fov = temp2
            r.ReadLine() ' water shift
            temp = r.ReadLine() ' angulation
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.angulation_midslice = temp2
            temp = r.ReadLine                       '.   offcenter                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.off_center_midslice = temp2
            r.ReadLine() ' flow
            r.ReadLine() ' presaturation
            r.ReadLine() ' phase encoding
            r.ReadLine() ' mtc
            r.ReadLine() ' spir
            r.ReadLine() ' epi
            r.ReadLine() ' ds
            r.ReadLine() ' diff
            r.ReadLine() ' dif echo time

            'now load info about the rec images
            r.ReadLine() '#
            r.ReadLine() '# === PIXEL VALUES =============================================================
            r.ReadLine() '#  PV = pixel value in REC file, FP = floating point value, DV = displayed value on console
            r.ReadLine() '#  RS = rescale slope,           RI = rescale intercept,    SS = scale slope
            r.ReadLine() '#  DV = PV * RS + RI             FP = DV / (RS * SS)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION DEFINITION =============================================
            r.ReadLine() '#  The rest of this file contains ONE line per image, this line contains the following information:
            r.ReadLine() '#
            r.ReadLine() '#  slice number      	                       (integer)
            r.ReadLine() '#  echo number                              (integer)
            r.ReadLine() '#  dynamic scan number                      (integer)
            r.ReadLine() '#  cardiac phase number                     (integer)
            r.ReadLine() '#  image_type_mr                            (integer)
            r.ReadLine() '#  scanning sequence                        (integer)
            r.ReadLine() '#  index in REC file (in images)            (integer)
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()

            r.ReadLine() '#  rescale intercept                        (float)
            r.ReadLine() '#  rescale slope                            (float)
            r.ReadLine() '#  scale slope                              (float)
            r.ReadLine() '#  window center                            (integer)
            r.ReadLine() '#  window width                             (integer)
            r.ReadLine() '#  image angulation (ap,fh,rl in degrees )  (3*float)
            r.ReadLine() '#'  image offcentre (ap,fh,rl in mm )        (3*float)

            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()

            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION ==========================================================
            r.ReadLine() '#sl ec dyn ph ty  idx (re)scale             window       angulation        offcentre         info     spacing   echo  dtime ttime diff  flip'
            r.ReadLine() '..

            'fill up the slices combo box and data in the array
            Dim i As Integer = 0
            Dim array_list As New ArrayList
            temp = r.ReadLine.ToString

            Dim start As Integer = 0
            Dim temp3 As String
            fov = New CoordinateD(0, 0, 0)

            Do While temp.Length > 100
                'recc_info(i) = New RECimagesInfo
                'declare array of objects and fill it
                Dim list_data As New RECimagesInfo
                list_data = New RECimagesInfo
                start = 0
                list_data.slice_number = CInt(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)

                list_data.image_type_mr = CInt(ReadNextAfterSpace(temp, start))
                list_data.index_in_rec_file = CInt(ReadNextAfterSpace(temp, start))

                par_struct.image_pixel_size = CInt(ReadNextAfterSpace(temp, start))
                pixel_size = CInt(par_struct.image_pixel_size)
                temp3 = ReadNextAfterSpace(temp, start)
                par_struct.recon_resolution_x_y(0) = CInt(ReadNextAfterSpace(temp, start))
                par_struct.recon_resolution_x_y(1) = CInt(ReadNextAfterSpace(temp, start))

                list_data.rescale_intercept = CSng(ReadNextAfterSpace(temp, start))
                list_data.rescale_slope = CSng(ReadNextAfterSpace(temp, start))
                list_data.scale_slope = CSng(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start) 'window center
                temp3 = ReadNextAfterSpace(temp, start) 'window width
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                par_struct.slice_thickness = CSng(ReadNextAfterSpace(temp, start)) 'slice thickness
                par_struct.slice_gap = CSng(ReadNextAfterSpace(temp, start)) 'slice gap
                temp3 = ReadNextAfterSpace(temp, start) 'display orientation
                list_data.slice_orientation = CInt(ReadNextAfterSpace(temp, start)) ' slice orientation
                temp3 = ReadNextAfterSpace(temp, start) 'fmri status indication
                temp3 = ReadNextAfterSpace(temp, start) 'image type

                fov.x = Convert.ToDouble(par_struct.recon_resolution_x_y(0) * Double.Parse(ReadNextAfterSpace(temp, start))) 'pixelspacingx
                fov.y = Convert.ToDouble(par_struct.recon_resolution_x_y(1) * Double.Parse(ReadNextAfterSpace(temp, start))) 'pixelspacingx

                array_list.Add(list_data)
                temp = r.ReadLine.ToString
                i += 1
            Loop
            par_struct.max_no_of_slices = i
            'now that we know the volume, move all the data to the recc_info structure and dispose the array_list
            ReDim recc_info(par_struct.max_no_of_slices)
            Dim a As Integer
            For a = 0 To par_struct.max_no_of_slices * 1 - 1
                recc_info(a) = CType(array_list.Item(a), RECimagesInfo)
            Next
            array_list.Clear()
            '            fov.x = CShort(CInt(par_struct.fov(2)))
            '            fov.y = CShort(CInt(par_struct.fov(0)))
            fov.z = (par_struct.slice_gap + par_struct.slice_thickness) * par_struct.max_no_of_slices
        End If
        'compute the actual fov
        rs = recc_info(0).rescale_slope
        ri = recc_info(0).rescale_intercept
        Dim center As New CoordinateD(0, 0, 0)
        Dim centerangle As New CoordinateD(0, 0, 0)

        
        center.x = par_struct.off_center_midslice(2)
        center.y = par_struct.off_center_midslice(0)
        center.z = par_struct.off_center_midslice(1)

        centerangle.x = Convert.ToDouble(par_struct.angulation_midslice(2))
        centerangle.y = Convert.ToDouble(par_struct.angulation_midslice(0))
        centerangle.z = Convert.ToDouble(par_struct.angulation_midslice(1))

        dimensions = New Coordinate(par_struct.recon_resolution_x_y(0), par_struct.recon_resolution_x_y(1), par_struct.max_no_of_slices)

        origin = center
        originAngle = centerangle
        orientation = recc_info(0).slice_orientation
        If orientation = 2 Then
            'sagittal slices
            Dim tt As Double = fov.x
            fov.x = fov.z
            fov.z = tt
        End If
        If orientation = 3 Then

        End If
    End Sub

    Public Function OpenOptimizedREC(ByVal fname As String, ByVal dimensions As Coordinate, ByVal pixel_size As Integer, ByVal rs As Single, ByVal ri As Single, ByVal orientation As Integer) As Matrix3DInt16
        'data() has X and Y interchanged
        Dim fn_img As String
        fn_img = Mid(fname, 1, Len(fname) - 3) & "rec"
        Dim fs As FileStream 'Declare new stream
        '    Try 'check if data.dat exsists.  If not error close
        fs = New FileStream(fn_img, FileMode.Open)
        Dim r As BinaryReader = New BinaryReader(fs) 'create reader
        Dim i As Long = 0
        Dim size As Int32 = (Convert.ToInt32(dimensions.x) * Convert.ToInt32(dimensions.y) * Convert.ToInt32(dimensions.z * 1))
        Dim data_list(size - 1, 0, 0) As Integer

        Select Case pixel_size
            Case 8
                Do While i < fs.Length 'reading one byte at a time
                    data_list(CInt(i), 0, 0) = CInt((Convert.ToInt16(r.ReadByte()) * rs) + ri)
                    i = i + 1
                Loop
            Case 16
                Do While i < (fs.Length / 2) 'because we are reading 2 bytes at a time
                    data_list(CInt(Convert.ToUInt64(i)), 0, 0) = CInt((r.ReadInt16 * rs) + ri)
                    i = i + 1
                Loop
        End Select
        'reading completed, REC files are always 4D matrices, although for single volume(Anatomy mostly), they can be converted to 3D matrices

        Dim data As New Matrix3DInt16(dimensions.x, dimensions.y, dimensions.z)
        data = reshape(data_list, size, 1, 1, dimensions.x, dimensions.y, dimensions.z)

        'reshape everything to Transverse
        Select Case orientation
            Case 2
                'sagittal
                Dim temp As New Matrix3DInt16(dimensions.z, dimensions.x, dimensions.y)
                For a As Integer = 0 To dimensions.z - 1
                    For b As Integer = 0 To dimensions.x - 1
                        For c As Integer = 0 To dimensions.y - 1
                            temp(dimensions.z - 1 - a, b, dimensions.y - 1 - c) = data(b, c, a)
                        Next
                    Next
                Next
                data = temp
            Case 3
                'coronal
                Dim temp As New Matrix3DInt16(dimensions.x, dimensions.z, dimensions.y)
                For a As Integer = 0 To dimensions.x - 1
                    For b As Integer = 0 To dimensions.z - 1
                        For c As Integer = 0 To dimensions.y - 1
                            temp(a, b, dimensions.y - 1 - c) = data(a, c, b)
                        Next
                    Next
                Next
                data = temp
        End Select

        Return data

    End Function
    Public Function reshape(ByVal matrix(,,) As Integer, ByVal matrix_x As Long, ByVal matrix_y As Long, ByVal matrix_z As Long, ByVal x As Integer, ByVal y As Integer, ByVal z As Integer) As Matrix3DInt16
        'RESHAPE(X,M,N) returns the M-by-N matrix whose elements
        'are taken columnwise from X.  An error results if X does
        'not have M*N elements.

        'it works like this:
        'it reads the first column of the first slice
        'then to the 2nd and then 3rd column of the first slice
        'then proceeds to the next slice
        'then to the next volume
        Dim data As New Matrix3DInt16(x, y, z)

        Dim array_list((CInt(matrix_x) * CInt(matrix_y) * CInt(matrix_z))) As Integer
        If matrix_x * matrix_y * matrix_z <> x * y * z Then
            MsgBox("Invalid call of function")
            Return Nothing
        Else
            'this part reads column wise in a array_list
            Dim a, b, c As Integer
            Dim d As Long
            d = 0
            For c = 0 To CInt(matrix_z - 1)
                For b = 0 To CInt(matrix_y - 1)
                    For a = 0 To CInt(matrix_x - 1)
                        array_list(CInt(d)) = matrix(a, b, c)
                        d += 1
                    Next
                Next
            Next
            'this part writes column wise from the array_list
            Dim counter As Long = 0
                For c = 0 To CInt(z - 1)
                For b = 0 To CInt(y - 1)
                    For a = 0 To CInt(x - 1)
                        data.data(a, b, c) = CShort(array_list(CInt(counter)))
                        counter = counter + 1
                    Next
                Next
            Next
            Return data
        End If
    End Function
    Public Function LoadSlicesFromVolume(ByVal data As Matrix4DInt16, ByVal volume_no As Integer) As Matrix3DInt16
        Dim ans As New Matrix3DInt16(data.i, data.j, data.k)
        Dim i, j, k As Integer
        For i = 0 To data.i - 1
            For j = 0 To data.j - 1
                For k = 0 To data.k - 1
                    ans.data(i, j, k) = data.data(i, j, k, volume_no)
                Next
            Next
        Next
        Return ans
    End Function
    Function ReadNextAfterSpace(ByVal Line As String, ByRef Start As Integer) As String
        'Start is 0 based
        Dim Ans As String
        Dim temp() As Char = Line.ToCharArray
        While temp(Start) = " "
            Start += 1
        End While
        'This is actual Start
        Dim st As Integer = Start
        While temp(Start) <> " " And Start < Line.Length - 1
            Start += 1
        End While
        Ans = Mid(Line, st + 1, Start - st + 1).Trim
        Return Ans
    End Function
    Sub Dispose()
        Me.LoadedRawData.Matrix.Dispose()
    End Sub
End Class

Public Class MultipleOptimizedParRecLoader
    Sub Dispose()
        Me.LoadedRawData.Matrix.Dispose()
    End Sub
    Dim LoadedRawData As MultiVolume
    Public Sub New(ByVal FileName As String)
        LoadedRawData = New MultiVolume
        'read the header list
        Dim dti_slice_info(0) As RECimagesInfo
        Dim no_of_volumes As Integer
        Dim pixel_size As Integer
        Dim rs, ri As Single 'needed to produce display 
        OpenOptimizedPAR(FileName, dti_slice_info, LoadedRawData.Center, LoadedRawData.OffCenterAngle, LoadedRawData.nx_ny_nz, LoadedRawData.FOV, no_of_volumes, pixel_size, rs, ri)

        'origin, dimensions and fov are already set
        'nx ny nz as dimensions

        'load the data into a specific volume
        LoadedRawData.Matrix = OpenOptimizedREC(FileName, LoadedRawData.nx_ny_nz, pixel_size, no_of_volumes, rs, ri)
        LoadedRawData.NumberOfVolumes = no_of_volumes
        'return the object
    End Sub
    Public Function GetData() As MultiVolume
        Return LoadedRawData
    End Function
    Sub OpenOptimizedPAR(ByVal fname As String, ByRef recc_info() As RECimagesInfo, ByRef origin As CoordinateD, ByRef originangle As CoordinateD, ByRef dimensions As Coordinate, ByRef fov As CoordinateD, ByRef no_of_volumes As Integer, ByRef pixel_size As Integer, ByRef rs As Single, ByRef ri As Single)
        Dim fs As FileStream 'Declare new stream

        fs = New IO.FileStream(fname.ToString, FileMode.Open)
        Dim r As StreamReader = New StreamReader(fs) 'create reader
        'read data into structure
        Dim par_struct As New PARfileInfo
        Dim temp As String
        Dim temp2() As Single
        r.ReadLine()                            '# === DATA DESCRIPTION FILE ======================================================
        r.ReadLine()                            '#
        r.ReadLine()                            '# CAUTION - Investigational device.
        r.ReadLine()                            '# Limited by Federal Law to investigational use.
        r.ReadLine()                            '#
        temp = r.ReadLine()                     '# Dataset name: e:\intera\TH_E004_3_1
        'par_struct.dataset_name = Mid(temp, 17, Len(temp) - 17 + 1)
        r.ReadLine()                            '#
        temp = r.ReadLine()                            '# CLINICAL TRYOUT             Research image export tool     V3
        If temp.Substring(61, 2) = "V3" Then
            'MsgBox("V3")
            r.ReadLine()                            '#
            r.ReadLine()                            '# === GENERAL INFORMATION ========================================================
            r.ReadLine()                            '#
            temp = r.ReadLine                       '.    Patient(Name)                      :   TH_E004_11292004()
            'par_struct.patient_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination(Name)                  :   Eureka()
            'par_struct.examination_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Protocol name                      :   T1W/3D/TFE
            'par_struct.protocol_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination date/time              :   2004.11.29 / 15:27:35
            'par_struct.exam_time_date = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Acquisition nr                     :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Reconstruction nr                  :   1
            'par_struct.reconstruction_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan Duration [sec]                :   113
            'par_struct.scan_duration_secs = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Max. number of cardiac phases      :   1
            'par_struct.max_no_of_cardic_phases = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of echoes              :   1
            'par_struct.max_no_of_echoes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of slices/locations    :   60
            par_struct.max_no_of_slices = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of dynamics            :   1
            'par_struct.max_no_of_dynamics = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of mixes               :   1
            'par_struct.max_no_of_mixes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Image pixel size [8 or 16 bits]    :   16
            par_struct.image_pixel_size = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            pixel_size = CInt(par_struct.image_pixel_size)
            temp = r.ReadLine                       '.    Technique                          :   T1TFE()
            'par_struct.technique = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan mode                          :   3D
            'par_struct.scanmode = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan resolution  (x, y)            :   240  150
            'temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 2)
            'par_struct.scan_resolution_x_y(0) = CInt(temp2(0))
            'par_struct.scan_resolution_x_y(1) = CInt(temp2(1))
            temp = r.ReadLine                       '.    Scan percentage                    :   100
            'par_struct.scan_percentage = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Recon resolution (x, y)            :   256  256
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 2)
            par_struct.recon_resolution_x_y(0) = CInt(temp2(0))
            par_struct.recon_resolution_x_y(1) = CInt(temp2(1))
            temp = r.ReadLine                       '.    Number of averages                 :   1
            'par_struct.no_of_averages = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Repetition time [msec]             :   9.31  
            'par_struct.repetition_time = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    FOV (ap,fh,rl) [mm]                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.fov = temp2
            temp = r.ReadLine                       '.    Slice thickness [mm]               :   2.00
            par_struct.slice_thickness = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Slice gap [mm]                     :   0.00
            par_struct.slice_gap = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Water Fat shift [pixels]           :   3.05
            'par_struct.water_fat_shift = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Angulation midslice(ap,fh,rl)[degr]:   0.00  0.00  0.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.angulation_midslice = temp2
            temp = r.ReadLine                       '.    Off Centre midslice(ap,fh,rl) [mm] :   -0.16  45.35  0.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.off_center_midslice = temp2
            temp = r.ReadLine                       '.    Flow compensation <0=no 1=yes> ?   :   0
            'par_struct.flow_compensation = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Presaturation     <0=no 1=yes> ?   :   1
            'par_struct.pre_saturation = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Cardiac frequency                  :   0
            'par_struct.cardiac_frequency = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Min. RR interval                   :   0
            'par_struct.min_rr_interval = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. RR interval                   :   0
            'par_struct.max_rr_interval = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Phase encoding velocity [cm/sec]   :   0.00  0.00  0.00
            'temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            'par_struct.phase_encoding_velocity = temp2
            temp = r.ReadLine                       '.    MTC               <0=no 1=yes> ?   :   0
            'par_struct.mtc = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    SPIR              <0=no 1=yes> ?   :   0
            'par_struct.spir = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    EPI factor        <0,1=no EPI>     :   1
            'par_struct.epi_factor = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    TURBO factor      <0=no turbo>     :   150
            'par_struct.turbo_factor = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Dynamic scan      <0=no 1=yes> ?   :   0
            'par_struct.dynamic_scan = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Diffusion         <0=no 1=yes> ?   :   0
            'par_struct.diffusion = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Diffusion echo time [msec]         :   0.00
            'par_struct.diffusion_echo_time = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Inversion delay [msec]             :   0.00
            'par_struct.inversion_delay = CSng(Mid(temp, 45, Len(temp) - 45 + 1))

            'now load info about the rec images

            r.ReadLine() '#
            r.ReadLine() '# === PIXEL VALUES =============================================================
            r.ReadLine() '#  PV = pixel value in REC file, FP = floating point value, DV = displayed value on console
            r.ReadLine() '#  RS = rescale slope,           RI = rescale intercept,    SS = scale slope
            r.ReadLine() '#  DV = PV * RS + RI             FP = DV / (RS * SS)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION DEFINITION =============================================
            r.ReadLine() '#  The rest of this file contains ONE line per image, this line contains the following information:
            r.ReadLine() '#
            r.ReadLine() '#  slice number      	                       (integer)
            r.ReadLine() '#  echo number                              (integer)
            r.ReadLine() '#  dynamic scan number                      (integer)
            r.ReadLine() '#  cardiac phase number                     (integer)
            r.ReadLine() '#  image_type_mr                            (integer)
            r.ReadLine() '#  scanning sequence                        (integer)
            r.ReadLine() '#  index in REC file (in images)            (integer)
            r.ReadLine() '#  rescale intercept                        (float)
            r.ReadLine() '#  rescale slope                            (float)
            r.ReadLine() '#  scale slope                              (float)
            r.ReadLine() '#  window center                            (integer)
            r.ReadLine() '#  window width                             (integer)
            r.ReadLine() '#  image angulation (ap,fh,rl in degrees )  (3*float)
            r.ReadLine() '#'  image offcentre (ap,fh,rl in mm )        (3*float)
            r.ReadLine() 'r.ReadLine()#'  image_display_orientation                (integer)
            r.ReadLine() '#  slice orientation ( TRA/SAG/COR )        (integer)
            r.ReadLine() '#  fmri_status_indication                   (integer)
            r.ReadLine() '#  image_type_ed_es  (end diast/end syst)   (integer)
            r.ReadLine() '#  pixel spacing (x,y) (in mm)              (2*float)
            r.ReadLine() '#  echo_time                                (float)
            r.ReadLine() '#  dyn_scan_begin_time                      (float)
            r.ReadLine() '#  trigger_time                             (float)
            r.ReadLine() '#  diffusion_b_factor                       (float)
            r.ReadLine() '#  image_flip_angle (in degrees)            (float)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION ==========================================================
            r.ReadLine() '#sl ec dyn ph ty  idx (re)scale             window       angulation        offcentre         info     spacing   echo  dtime ttime diff  flip'
            r.ReadLine() '..

            'fill up the slices combo box and data in the array
            Dim i As Integer = 0
            Dim array_list As New ArrayList
            temp = r.ReadLine.ToString

            Dim start As Integer = 0
            Dim temp3 As String

            Do While temp.Length > 100
                'recc_info(i) = New RECimagesInfo
                'declare array of objects and fill it
                Dim list_data As New RECimagesInfo
                list_data = New RECimagesInfo
                start = 0
                list_data.slice_number = CInt(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)

                list_data.image_type_mr = CInt(ReadNextAfterSpace(temp, start))
                '  MsgBox(list_data.image_type_mr)
                list_data.index_in_rec_file = CInt(ReadNextAfterSpace(temp, start))
                list_data.rescale_intercept = CSng(ReadNextAfterSpace(temp, start))
                list_data.rescale_slope = CSng(ReadNextAfterSpace(temp, start))
                list_data.scale_slope = CSng(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                list_data.slice_orientation = CInt(ReadNextAfterSpace(temp, start))
                i = CInt(list_data.index_in_rec_file)

                array_list.Add(list_data)
                temp = r.ReadLine.ToString
            Loop
            no_of_volumes = CInt((i + 1) / par_struct.max_no_of_slices)

            'now that we know the volume, move all the data to the recc_info structure and dispose the array_list
            ReDim recc_info(par_struct.max_no_of_slices * no_of_volumes)
            Dim a As Integer
            For a = 0 To par_struct.max_no_of_slices * no_of_volumes - 1
                recc_info(a) = CType(array_list.Item(a), RECimagesInfo)
            Next
            array_list.Clear()
        ElseIf temp.Substring(61, 2) = "V4" Then
            'MsgBox("V4")
            r.ReadLine() '#
            r.ReadLine()                            '# === GENERAL INFORMATION ========================================================
            r.ReadLine()                            '#
            temp = r.ReadLine                       '.    Patient(Name)                      :   TH_E004_11292004()
            'par_struct.patient_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination(Name)                  :   Eureka()
            'par_struct.examination_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Protocol name                      :   T1W/3D/TFE
            'par_struct.protocol_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination date/time              :   2004.11.29 / 15:27:35
            'par_struct.exam_time_date = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Series Type                      :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Acquisition nr                     :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Reconstruction nr                  :   1
            'par_struct.reconstruction_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan Duration [sec]                :   113
            'par_struct.scan_duration_secs = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Max. number of cardiac phases      :   1
            'par_struct.max_no_of_cardic_phases = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of echoes              :   1
            'par_struct.max_no_of_echoes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of slices/locations    :   60
            par_struct.max_no_of_slices = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of dynamics            :   1
            'par_struct.max_no_of_dynamics = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of mixes               :   1
            'par_struct.max_no_of_mixes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            r.ReadLine() 'Patient Position
            r.ReadLine() 'Preparation Direction
            r.ReadLine() ' technique
            r.ReadLine() ' scan resolution
            r.ReadLine() 'scan mode
            r.ReadLine() 'repitition time
            temp = r.ReadLine                       '.    FOV (ap,fh,rl) [mm]                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.fov = temp2
            r.ReadLine() ' water shift
            temp = r.ReadLine() ' angulation
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.angulation_midslice = temp2
            temp = r.ReadLine                       '.   offcenter                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.off_center_midslice = temp2
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()

            'now load info about the rec images

            r.ReadLine() '#
            r.ReadLine() '# === PIXEL VALUES =============================================================
            r.ReadLine() '#  PV = pixel value in REC file, FP = floating point value, DV = displayed value on console
            r.ReadLine() '#  RS = rescale slope,           RI = rescale intercept,    SS = scale slope
            r.ReadLine() '#  DV = PV * RS + RI             FP = DV / (RS * SS)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION DEFINITION =============================================
            r.ReadLine() '#  The rest of this file contains ONE line per image, this line contains the following information:
            r.ReadLine() '#
            r.ReadLine() '#  slice number      	                       (integer)
            r.ReadLine() '#  echo number                              (integer)
            r.ReadLine() '#  dynamic scan number                      (integer)
            r.ReadLine() '#  cardiac phase number                     (integer)
            r.ReadLine() '#  image_type_mr                            (integer)
            r.ReadLine() '#  scanning sequence                        (integer)
            r.ReadLine() '#  index in REC file (in images)            (integer)
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()

            r.ReadLine() '#  rescale intercept                        (float)
            r.ReadLine() '#  rescale slope                            (float)
            r.ReadLine() '#  scale slope                              (float)
            r.ReadLine() '#  window center                            (integer)
            r.ReadLine() '#  window width                             (integer)
            r.ReadLine() '#  image angulation (ap,fh,rl in degrees )  (3*float)
            r.ReadLine() '#'  image offcentre (ap,fh,rl in mm )        (3*float)

            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()

            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION ==========================================================
            r.ReadLine() '#sl ec dyn ph ty  idx (re)scale             window       angulation        offcentre         info     spacing   echo  dtime ttime diff  flip'
            r.ReadLine() '..

            'fill up the slices combo box and data in the array
            Dim i As Integer = 0
            Dim array_list As New ArrayList
            temp = r.ReadLine.ToString

            Dim start As Integer = 0
            Dim temp3 As String
            Do While temp.Length > 100
                'recc_info(i) = New RECimagesInfo
                'declare array of objects and fill it
                Dim list_data As New RECimagesInfo
                list_data = New RECimagesInfo
                start = 0
                list_data.slice_number = CInt(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)

                list_data.image_type_mr = CInt(ReadNextAfterSpace(temp, start))
                list_data.index_in_rec_file = CInt(ReadNextAfterSpace(temp, start))

                par_struct.image_pixel_size = CInt(ReadNextAfterSpace(temp, start))
                pixel_size = CInt(par_struct.image_pixel_size)
                temp3 = ReadNextAfterSpace(temp, start)
                par_struct.recon_resolution_x_y(0) = CInt(ReadNextAfterSpace(temp, start))
                par_struct.recon_resolution_x_y(1) = CInt(ReadNextAfterSpace(temp, start))
                list_data.rescale_intercept = CSng(ReadNextAfterSpace(temp, start))
                list_data.rescale_slope = CSng(ReadNextAfterSpace(temp, start))
                list_data.scale_slope = CSng(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start) 'window center
                temp3 = ReadNextAfterSpace(temp, start) 'window width
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                par_struct.slice_thickness = CSng(ReadNextAfterSpace(temp, start)) 'slice thickness
                par_struct.slice_gap = CSng(ReadNextAfterSpace(temp, start)) 'slice gap
                temp3 = ReadNextAfterSpace(temp, start) 'display orientation
                list_data.slice_orientation = CInt(ReadNextAfterSpace(temp, start)) ' slice orientation

                i = CInt(list_data.index_in_rec_file)
                array_list.Add(list_data)
                temp = r.ReadLine.ToString
            Loop
            no_of_volumes = CInt((i + 1) / par_struct.max_no_of_slices)

            'now that we know the volume, move all the data to the recc_info structure and dispose the array_list
            ReDim recc_info(par_struct.max_no_of_slices * no_of_volumes)
            Dim a As Integer
            For a = 0 To par_struct.max_no_of_slices * no_of_volumes - 1
                recc_info(a) = CType(array_list.Item(a), RECimagesInfo)
            Next
            array_list.Clear()
        End If

        'compute the actual fov
        rs = recc_info(0).rescale_slope
        ri = recc_info(0).rescale_intercept
        fov = New CoordinateD(0, 0, 0)

        Dim center As New CoordinateD(0, 0, 0)
        Dim centerangle As New CoordinateD(0, 0, 0)

        fov.x = (par_struct.fov(2))
        fov.y = (par_struct.fov(0))
        fov.z = (par_struct.fov(1)) '(par_struct.slice_gap + par_struct.slice_thickness) * par_struct.max_no_of_slices

        center.x = (par_struct.off_center_midslice(2))
        center.y = (par_struct.off_center_midslice(0))
        center.z = (par_struct.off_center_midslice(1))

        centerangle.x = (par_struct.angulation_midslice(2))
        centerangle.y = (par_struct.angulation_midslice(0))
        centerangle.z = (par_struct.angulation_midslice(1))

        Select Case recc_info(0).slice_orientation
            Case 1
                'slice orientation transversal
                'xy seen
            Case 2
                'yz seen
                MsgBox("Sagittal Slice orientation in Diffusion" & _
                vbNewLine & "Weighted data is not supported.")
            Case 3
                MsgBox("Coronal Slice orientation in Diffusion" & _
                vbNewLine & "Weighted data is not supported.")
        End Select
        'MsgBox("The recommended resize for FA is : " & vbNewLine & _
        '        IIf(recc_info(0).slice_orientation = 1, "Transverse slice width ratio: " & vbTab & vbTab & fov.x _
        '        & vbNewLine & "Transverse slice height ratio: " & vbTab & vbTab & fov.y _
        '        & vbNewLine & "Coronal slice height ratio: " & vbTab & vbTab & fov.z _
        '        & vbNewLine & "Center coordinate point: " & vbTab & vbTab & center.x & "," & center.y & "," & center.z _
        '        , "") & _
        '                ".", MsgBoxStyle.Information, "Slice Information")

        origin = center
        originangle = centerangle
        dimensions = New Coordinate(par_struct.recon_resolution_x_y(0), par_struct.recon_resolution_x_y(1), par_struct.max_no_of_slices)
    End Sub
    Function ReadNextAfterSpace(ByVal Line As String, ByRef Start As Integer) As String
        'Start is 0 based
        Dim Ans As String
        Dim temp() As Char = Line.ToCharArray
        While temp(Start) = " "
            Start += 1
        End While
        'This is actual Start
        Dim st As Integer = Start
        While temp(Start) <> " " And Start < Line.Length - 1
            Start += 1
        End While
        Ans = Mid(Line, st + 1, Start - st + 1).Trim
        Return Ans
    End Function
    Public Function OpenOptimizedREC(ByVal fname As String, ByVal dimensions As Coordinate, ByVal pixel_size As Integer, ByVal volumes As Integer, ByVal rs As Single, ByVal ri As Single) As Matrix4DInt16
        'data() has X and Y interchanged
        Dim fn_img As String
        fn_img = Mid(fname, 1, Len(fname) - 3) & "rec"
        Dim fs As FileStream 'Declare new stream
        '    Try 'check if data.dat exsists.  If not error close
        fs = New FileStream(fn_img, FileMode.Open)
        Dim r As BinaryReader = New BinaryReader(fs) 'create reader
        Dim i As Int32 = 0
        Dim size As Int32 = (Convert.ToInt32(dimensions.x) * Convert.ToInt32(dimensions.y) * Convert.ToInt32(dimensions.z * volumes))
        Dim data_list(size - 1, 0, 0) As Integer

        Select Case pixel_size
            Case 8
                Do While i < fs.Length 'reading one byte at a time
                    data_list(CInt(i), 0, 0) = CInt((Convert.ToInt16(r.ReadByte()) * rs) + ri)
                    i = i + 1
                Loop
            Case 16
                Do While i < (fs.Length / 2) 'because we are reading 2 bytes at a time
                    'MsgBox(r.PeekChar)
                    data_list(CInt(i), 0, 0) = CInt((r.ReadInt16)) ' * rs) + ri)
                    If data_list(CInt(i), 0, 0) < 0 Then MsgBox("caught")
                    i = i + 1
                Loop
        End Select
        'reading completed, REC files are always 4D matrices, although for single volume(Anatomy mostly), they can be converted to 3D matrices

        Dim data As New Matrix4DInt16(dimensions.y, dimensions.x, dimensions.z, volumes)
        data = reshape(data_list, size, 1, 1, dimensions.y, dimensions.x, dimensions.z, volumes)
        Return data

    End Function
    Public Function LoadSlicesFromVolume(ByVal data As Matrix4DInt16, ByVal volume_no As Integer) As Matrix3DInt16
        Dim ans As New Matrix3DInt16(data.i, data.j, data.k)
        Dim i, j, k As Integer
        For i = 0 To data.i - 1
            For j = 0 To data.j - 1
                For k = 0 To data.k - 1
                    ans.data(i, j, k) = data.data(i, j, k, volume_no)
                Next
            Next
        Next
        Return ans
    End Function
End Class

Public Class PerfusionParRecLoader
    Sub Dispose()
        Me.LoadedRawData.Matrix.Dispose()
    End Sub
    Dim LoadedRawData As MultiVolume
    Dim IndexArray() As Integer

    Public Sub New(ByVal FileName As String)
        LoadedRawData = New MultiVolume
        'read the header list
        Dim dti_slice_info(0) As RECimagesInfo
        Dim no_of_volumes As Integer
        Dim pixel_size As Integer
        Dim rs, ri As Single 'needed to produce display 
        OpenOptimizedPAR(FileName, dti_slice_info, LoadedRawData.Center, LoadedRawData.OffCenterAngle, LoadedRawData.nx_ny_nz, LoadedRawData.FOV, no_of_volumes, pixel_size, rs, ri, IndexArray)

        'origin, dimensions and fov are already set
        'nx ny nz as dimensions

        'load the data into a specific volume
        LoadedRawData.Matrix = OpenOptimizedREC(FileName, LoadedRawData.nx_ny_nz, pixel_size, no_of_volumes, rs, ri)
        LoadedRawData.NumberOfVolumes = no_of_volumes
        'return the object
    End Sub
    Public Function GetData() As MultiVolume
        Return LoadedRawData
    End Function
    Sub OpenOptimizedPAR(ByVal fname As String, ByRef recc_info() As RECimagesInfo, ByRef origin As CoordinateD, ByRef originAngle As CoordinateD, ByRef dimensions As Coordinate, ByRef fov As CoordinateD, ByRef no_of_volumes As Integer, ByRef pixel_size As Integer, ByRef rs As Single, ByRef ri As Single, ByRef indexarray() As Integer)
        Dim fs As FileStream 'Declare new stream

        fs = New IO.FileStream(fname.ToString, FileMode.Open)
        Dim r As StreamReader = New StreamReader(fs) 'create reader
        'read data into structure
        Dim par_struct As New PARfileInfo
        Dim temp As String
        Dim temp2() As Single
        r.ReadLine()                            '# === DATA DESCRIPTION FILE ======================================================
        r.ReadLine()                            '#
        r.ReadLine()                            '# CAUTION - Investigational device.
        r.ReadLine()                            '# Limited by Federal Law to investigational use.
        r.ReadLine()                            '#
        temp = r.ReadLine()                     '# Dataset name: e:\intera\TH_E004_3_1
        'par_struct.dataset_name = Mid(temp, 17, Len(temp) - 17 + 1)
        r.ReadLine()   '#
        temp = r.ReadLine()                            '# CLINICAL TRYOUT             Research image export tool     V3
        If temp.Contains("V3") Then
            MsgBox("V3")
            r.ReadLine()                            '#
            r.ReadLine()                            '# === GENERAL INFORMATION ========================================================
            r.ReadLine()                            '#
            temp = r.ReadLine                       '.    Patient(Name)                      :   TH_E004_11292004()
            'par_struct.patient_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination(Name)                  :   Eureka()
            'par_struct.examination_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Protocol name                      :   T1W/3D/TFE
            'par_struct.protocol_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination date/time              :   2004.11.29 / 15:27:35
            'par_struct.exam_time_date = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Acquisition nr                     :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Reconstruction nr                  :   1
            'par_struct.reconstruction_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan Duration [sec]                :   113
            'par_struct.scan_duration_secs = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Max. number of cardiac phases      :   1
            'par_struct.max_no_of_cardic_phases = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of echoes              :   1
            'par_struct.max_no_of_echoes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of slices/locations    :   60
            par_struct.max_no_of_slices = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of dynamics            :   1
            'par_struct.max_no_of_dynamics = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of mixes               :   1
            'par_struct.max_no_of_mixes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Image pixel size [8 or 16 bits]    :   16
            par_struct.image_pixel_size = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            pixel_size = CInt(par_struct.image_pixel_size)
            temp = r.ReadLine                       '.    Technique                          :   T1TFE()
            'par_struct.technique = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan mode                          :   3D
            'par_struct.scanmode = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan resolution  (x, y)            :   240  150
            'temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 2)
            'par_struct.scan_resolution_x_y(0) = CInt(temp2(0))
            'par_struct.scan_resolution_x_y(1) = CInt(temp2(1))
            temp = r.ReadLine                       '.    Scan percentage                    :   100
            'par_struct.scan_percentage = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Recon resolution (x, y)            :   256  256
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 2)
            par_struct.recon_resolution_x_y(0) = CInt(temp2(0))
            par_struct.recon_resolution_x_y(1) = CInt(temp2(1))
            temp = r.ReadLine                       '.    Number of averages                 :   1
            'par_struct.no_of_averages = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Repetition time [msec]             :   9.31  
            'par_struct.repetition_time = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    FOV (ap,fh,rl) [mm]                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.fov = temp2
            temp = r.ReadLine                       '.    Slice thickness [mm]               :   2.00
            par_struct.slice_thickness = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Slice gap [mm]                     :   0.00
            par_struct.slice_gap = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Water Fat shift [pixels]           :   3.05
            'par_struct.water_fat_shift = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Angulation midslice(ap,fh,rl)[degr]:   0.00  0.00  0.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.angulation_midslice = temp2
            temp = r.ReadLine                       '.    Off Centre midslice(ap,fh,rl) [mm] :   -0.16  45.35  0.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.off_center_midslice = temp2
            temp = r.ReadLine                       '.    Flow compensation <0=no 1=yes> ?   :   0
            'par_struct.flow_compensation = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Presaturation     <0=no 1=yes> ?   :   1
            'par_struct.pre_saturation = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Cardiac frequency                  :   0
            'par_struct.cardiac_frequency = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Min. RR interval                   :   0
            'par_struct.min_rr_interval = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. RR interval                   :   0
            'par_struct.max_rr_interval = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Phase encoding velocity [cm/sec]   :   0.00  0.00  0.00
            'temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            'par_struct.phase_encoding_velocity = temp2
            temp = r.ReadLine                       '.    MTC               <0=no 1=yes> ?   :   0
            'par_struct.mtc = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    SPIR              <0=no 1=yes> ?   :   0
            'par_struct.spir = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    EPI factor        <0,1=no EPI>     :   1
            'par_struct.epi_factor = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    TURBO factor      <0=no turbo>     :   150
            'par_struct.turbo_factor = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Dynamic scan      <0=no 1=yes> ?   :   0
            'par_struct.dynamic_scan = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Diffusion         <0=no 1=yes> ?   :   0
            'par_struct.diffusion = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Diffusion echo time [msec]         :   0.00
            'par_struct.diffusion_echo_time = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Inversion delay [msec]             :   0.00
            'par_struct.inversion_delay = CSng(Mid(temp, 45, Len(temp) - 45 + 1))

            'now load info about the rec images

            r.ReadLine() '#
            r.ReadLine() '# === PIXEL VALUES =============================================================
            r.ReadLine() '#  PV = pixel value in REC file, FP = floating point value, DV = displayed value on console
            r.ReadLine() '#  RS = rescale slope,           RI = rescale intercept,    SS = scale slope
            r.ReadLine() '#  DV = PV * RS + RI             FP = DV / (RS * SS)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION DEFINITION =============================================
            r.ReadLine() '#  The rest of this file contains ONE line per image, this line contains the following information:
            r.ReadLine() '#
            r.ReadLine() '#  slice number      	                       (integer)
            r.ReadLine() '#  echo number                              (integer)
            r.ReadLine() '#  dynamic scan number                      (integer)
            r.ReadLine() '#  cardiac phase number                     (integer)
            r.ReadLine() '#  image_type_mr                            (integer)
            r.ReadLine() '#  scanning sequence                        (integer)
            r.ReadLine() '#  index in REC file (in images)            (integer)
            r.ReadLine() '#  rescale intercept                        (float)
            r.ReadLine() '#  rescale slope                            (float)
            r.ReadLine() '#  scale slope                              (float)
            r.ReadLine() '#  window center                            (integer)
            r.ReadLine() '#  window width                             (integer)
            r.ReadLine() '#  image angulation (ap,fh,rl in degrees )  (3*float)
            r.ReadLine() '#'  image offcentre (ap,fh,rl in mm )        (3*float)
            r.ReadLine() 'r.ReadLine()#'  image_display_orientation                (integer)
            r.ReadLine() '#  slice orientation ( TRA/SAG/COR )        (integer)
            r.ReadLine() '#  fmri_status_indication                   (integer)
            r.ReadLine() '#  image_type_ed_es  (end diast/end syst)   (integer)
            r.ReadLine() '#  pixel spacing (x,y) (in mm)              (2*float)
            r.ReadLine() '#  echo_time                                (float)
            r.ReadLine() '#  dyn_scan_begin_time                      (float)
            r.ReadLine() '#  trigger_time                             (float)
            r.ReadLine() '#  diffusion_b_factor                       (float)
            r.ReadLine() '#  image_flip_angle (in degrees)            (float)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION ==========================================================
            r.ReadLine() '#sl ec dyn ph ty  idx (re)scale             window       angulation        offcentre         info     spacing   echo  dtime ttime diff  flip'
            r.ReadLine() '..

            'fill up the slices combo box and data in the array
            Dim i As Integer = 0
            Dim array_list As New ArrayList
            temp = r.ReadLine.ToString

            Dim start As Integer = 0
            Dim temp3 As String

            Do While temp.Length > 100
                'recc_info(i) = New RECimagesInfo
                'declare array of objects and fill it
                Dim list_data As New RECimagesInfo
                list_data = New RECimagesInfo
                start = 0
                list_data.slice_number = CInt(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)

                list_data.image_type_mr = CInt(ReadNextAfterSpace(temp, start))
                list_data.index_in_rec_file = CInt(ReadNextAfterSpace(temp, start))
                list_data.rescale_intercept = CSng(ReadNextAfterSpace(temp, start))
                list_data.rescale_slope = CSng(ReadNextAfterSpace(temp, start))
                list_data.scale_slope = CSng(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                list_data.slice_orientation = CInt(ReadNextAfterSpace(temp, start))
                i = CInt(list_data.index_in_rec_file)
                ReDim Preserve indexarray(i)
                indexarray(i) = list_data.image_type_mr

                array_list.Add(list_data)
                temp = r.ReadLine.ToString
            Loop
            no_of_volumes = CInt((i + 1) / (par_struct.max_no_of_slices * 2))

            'now that we know the volume, move all the data to the recc_info structure and dispose the array_list
            ReDim recc_info(par_struct.max_no_of_slices * no_of_volumes)
            Dim a As Integer
            For a = 0 To par_struct.max_no_of_slices * no_of_volumes - 1
                recc_info(a) = CType(array_list.Item(a), RECimagesInfo)
            Next
            array_list.Clear()
        ElseIf temp.Contains("V4") Then
            r.ReadLine() '#
            r.ReadLine()                            '# === GENERAL INFORMATION ========================================================
            r.ReadLine()                            '#
            temp = r.ReadLine                       '.    Patient(Name)                      :   TH_E004_11292004()
            'par_struct.patient_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination(Name)                  :   Eureka()
            'par_struct.examination_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Protocol name                      :   T1W/3D/TFE
            'par_struct.protocol_name = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Examination date/time              :   2004.11.29 / 15:27:35
            'par_struct.exam_time_date = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Series Type                      :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Acquisition nr                     :   3
            'par_struct.acquisition_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Reconstruction nr                  :   1
            'par_struct.reconstruction_nr = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Scan Duration [sec]                :   113
            'par_struct.scan_duration_secs = Mid(temp, 45, Len(temp) - 45 + 1)
            temp = r.ReadLine                       '.    Max. number of cardiac phases      :   1
            'par_struct.max_no_of_cardic_phases = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of echoes              :   1
            'par_struct.max_no_of_echoes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of slices/locations    :   60
            par_struct.max_no_of_slices = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of dynamics            :   1
            'par_struct.max_no_of_dynamics = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            temp = r.ReadLine                       '.    Max. number of mixes               :   1
            'par_struct.max_no_of_mixes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
            r.ReadLine() 'Patient Position
            r.ReadLine() 'Preparation Direction
            r.ReadLine() ' technique
            r.ReadLine() ' scan resolution
            r.ReadLine() 'scan mode
            r.ReadLine() 'repitition time
            temp = r.ReadLine                       '.    FOV (ap,fh,rl) [mm]                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.fov = temp2
            r.ReadLine() ' water shift
            temp = r.ReadLine() ' angulation
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.angulation_midslice = temp2
            temp = r.ReadLine                       '.   offcenter                :   230.00  120.00  230.00
            temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
            par_struct.off_center_midslice = temp2
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()

            'now load info about the rec images
            r.ReadLine() '#
            r.ReadLine() '# === PIXEL VALUES =============================================================
            r.ReadLine() '#  PV = pixel value in REC file, FP = floating point value, DV = displayed value on console
            r.ReadLine() '#  RS = rescale slope,           RI = rescale intercept,    SS = scale slope
            r.ReadLine() '#  DV = PV * RS + RI             FP = DV / (RS * SS)
            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION DEFINITION =============================================
            r.ReadLine() '#  The rest of this file contains ONE line per image, this line contains the following information:
            r.ReadLine() '#
            r.ReadLine() '#  slice number      	                       (integer)
            r.ReadLine() '#  echo number                              (integer)
            r.ReadLine() '#  dynamic scan number                      (integer)
            r.ReadLine() '#  cardiac phase number                     (integer)
            r.ReadLine() '#  image_type_mr                            (integer)
            r.ReadLine() '#  scanning sequence                        (integer)
            r.ReadLine() '#  index in REC file (in images)            (integer)
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()

            r.ReadLine() '#  rescale intercept                        (float)
            r.ReadLine() '#  rescale slope                            (float)
            r.ReadLine() '#  scale slope                              (float)
            r.ReadLine() '#  window center                            (integer)
            r.ReadLine() '#  window width                             (integer)
            r.ReadLine() '#  image angulation (ap,fh,rl in degrees )  (3*float)
            r.ReadLine() '#'  image offcentre (ap,fh,rl in mm )        (3*float)

            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()
            r.ReadLine()

            r.ReadLine() '#
            r.ReadLine() '# === IMAGE INFORMATION ==========================================================
            r.ReadLine() '#sl ec dyn ph ty  idx (re)scale             window       angulation        offcentre         info     spacing   echo  dtime ttime diff  flip'
            r.ReadLine() '..

            'fill up the slices combo box and data in the array
            Dim i As Integer = 0
            Dim array_list As New ArrayList
            temp = r.ReadLine.ToString

            Dim start As Integer = 0
            Dim temp3 As String
            Do While temp.Length > 100
                'recc_info(i) = New RECimagesInfo
                'declare array of objects and fill it
                Dim list_data As New RECimagesInfo
                list_data = New RECimagesInfo
                start = 0
                list_data.slice_number = CInt(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)
                temp3 = ReadNextAfterSpace(temp, start)

                list_data.image_type_mr = CInt(ReadNextAfterSpace(temp, start))
                list_data.index_in_rec_file = CInt(ReadNextAfterSpace(temp, start))

                par_struct.image_pixel_size = CInt(ReadNextAfterSpace(temp, start))
                pixel_size = CInt(par_struct.image_pixel_size)
                temp3 = ReadNextAfterSpace(temp, start)
                par_struct.recon_resolution_x_y(0) = CInt(ReadNextAfterSpace(temp, start))
                par_struct.recon_resolution_x_y(1) = CInt(ReadNextAfterSpace(temp, start))
                list_data.rescale_intercept = CSng(ReadNextAfterSpace(temp, start))
                list_data.rescale_slope = CSng(ReadNextAfterSpace(temp, start))
                list_data.scale_slope = CSng(ReadNextAfterSpace(temp, start))
                temp3 = ReadNextAfterSpace(temp, start) 'window center
                temp3 = ReadNextAfterSpace(temp, start) 'window width
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image angulation
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                temp3 = ReadNextAfterSpace(temp, start) 'image offcenter
                par_struct.slice_thickness = CSng(ReadNextAfterSpace(temp, start)) 'slice thickness
                par_struct.slice_gap = CSng(ReadNextAfterSpace(temp, start)) 'slice gap
                temp3 = ReadNextAfterSpace(temp, start) 'display orientation
                list_data.slice_orientation = CInt(ReadNextAfterSpace(temp, start)) ' slice orientation

                i = CInt(list_data.index_in_rec_file)
                ReDim Preserve indexarray(i)
                indexarray(i) = list_data.image_type_mr

                array_list.Add(list_data)
                temp = r.ReadLine.ToString

            Loop
            no_of_volumes = CInt((i + 1) / (par_struct.max_no_of_slices * 2))
            'now that we know the volume, move all the data to the recc_info structure and dispose the array_list
            ReDim recc_info(par_struct.max_no_of_slices * no_of_volumes)
            Dim a As Integer
            For a = 0 To par_struct.max_no_of_slices * no_of_volumes - 1
                recc_info(a) = CType(array_list.Item(a), RECimagesInfo)
            Next
            array_list.Clear()
        End If

        'compute the actual fov
        rs = recc_info(0).rescale_slope
        ri = recc_info(0).rescale_intercept
        fov = New CoordinateD(0, 0, 0)

        Dim center As New CoordinateD(0, 0, 0)
        Dim centerangle As New CoordinateD(0, 0, 0)

        fov.x = CShort(CInt(par_struct.fov(2)))
        fov.y = CShort(CInt(par_struct.fov(0)))
        fov.z = CShort(CInt(par_struct.fov(1))) '(par_struct.slice_gap + par_struct.slice_thickness) * par_struct.max_no_of_slices

        center.x = par_struct.off_center_midslice(2)
        center.y = par_struct.off_center_midslice(0)
        center.z = par_struct.off_center_midslice(1)

        centerangle.x = CShort(CInt(par_struct.angulation_midslice(2)))
        centerangle.y = CShort(CInt(par_struct.angulation_midslice(0)))
        centerangle.z = CShort(CInt(par_struct.angulation_midslice(1)))

        Select Case recc_info(0).slice_orientation
            Case 1
                'slice orientation transversal
                'xy seen
            Case 2
                'yz seen
                MsgBox("Sagittal Slice orientation in Perfusion" & _
                vbNewLine & "Weighted data is not supported.")
            Case 3
                MsgBox("Coronal Slice orientation in Perfusion" & _
                vbNewLine & "Weighted data is not supported.")
        End Select

        origin = center
        dimensions = New Coordinate(par_struct.recon_resolution_x_y(0), par_struct.recon_resolution_x_y(1), par_struct.max_no_of_slices * 2)
        originAngle = centerangle
    End Sub
    Function ReadNextAfterSpace(ByVal Line As String, ByRef Start As Integer) As String
        'Start is 0 based
        Dim Ans As String
        Dim temp() As Char = Line.ToCharArray
        While temp(Start) = " "
            Start += 1
        End While
        'This is actual Start
        Dim st As Integer = Start
        While temp(Start) <> " " And Start < Line.Length - 1
            Start += 1
        End While
        Ans = Mid(Line, st + 1, Start - st + 1).Trim
        Return Ans
    End Function
    Public Function OpenOptimizedREC(ByVal fname As String, ByVal dimensions As Coordinate, ByVal pixel_size As Integer, ByVal volumes As Integer, ByVal rs As Single, ByVal ri As Single) As Matrix4DInt16
        'data() has X and Y interchanged
        Dim fn_img As String
        fn_img = Mid(fname, 1, Len(fname) - 3) & "rec"
        Dim fs As FileStream 'Declare new stream
        '    Try 'check if data.dat exsists.  If not error close
        fs = New FileStream(fn_img, FileMode.Open)
        Dim r As BinaryReader = New BinaryReader(fs) 'create reader
        Dim i As Int32 = 0
        Dim size As Int32 = (Convert.ToInt32(dimensions.x) * Convert.ToInt32(dimensions.y) * Convert.ToInt32(dimensions.z * volumes))
        Dim data_list(size - 1, 0, 0) As Integer

        Select Case pixel_size
            Case 8
                Do While i < fs.Length 'reading one byte at a time
                    data_list(CInt(i), 0, 0) = CInt((Convert.ToInt16(r.ReadByte()) * rs) + ri)
                    i = i + 1
                Loop
            Case 16
                Do While i < (fs.Length / 2) 'because we are reading 2 bytes at a time
                    'MsgBox(r.PeekChar)
                    data_list(CInt(i), 0, 0) = CInt((r.ReadInt16)) ' * rs) + ri)
                    If data_list(CInt(i), 0, 0) < 0 Then MsgBox("caught")
                    i = i + 1
                Loop
        End Select
        'reading completed, REC files are always 4D matrices, although for single volume(Anatomy mostly), they can be converted to 3D matrices

        Dim data As New Matrix4DInt16(dimensions.y, dimensions.x, dimensions.z, volumes)
        data = reshape(data_list, size, 1, 1, dimensions.y, dimensions.x, dimensions.z, volumes)

        Return data

    End Function
    Public Function LoadSlicesFromVolume(ByVal data As Matrix4DInt16, ByVal volume_no As Integer) As Matrix3DInt16
        Dim ans As New Matrix3DInt16(data.i, data.j, data.k)
        Dim i, j, k As Integer
        For i = 0 To data.i - 1
            For j = 0 To data.j - 1
                For k = 0 To data.k - 1
                    ans.data(i, j, k) = data.data(i, j, k, volume_no)
                Next
            Next
        Next
        Return ans
    End Function
    'Public Function PartitionData() As Matrix3DSingle

    '    Dim Temp5D As New Matrix4Ddouble(Me.LoadedRawData.Matrix.i, Me.LoadedRawData.Matrix.j, Me.LoadedRawData.Matrix.k, 2)

    '    Dim i, j, k, l, counter As Integer
    '    Dim ctr1, ctr2 As Integer

    '    counter = 0

    '    For l = 0 To Me.LoadedRawData.Matrix.l - 1
    '        If l Mod 2 = 0 Then
    '            ctr1 = 0
    '            ctr2 = 0
    '        End If

    '        For k = 0 To Me.LoadedRawData.Matrix.k - 1
    '            'we are at slice level now
    '            If IndexArray(counter) = "1" Then
    '                'load into 1st multivolume
    '                For j = 0 To Me.LoadedRawData.Matrix.j - 1
    '                    For i = 0 To Me.LoadedRawData.Matrix.i - 1
    '                        Temp5D.data(i, j, ctr1, 0) += Me.LoadedRawData.Matrix.data(i, j, k, l)
    '                    Next
    '                Next
    '                ctr1 += 1

    '            ElseIf IndexArray(counter) = "4" Then
    '                'load into 2nd multivolume
    '                For j = 0 To Me.LoadedRawData.Matrix.j - 1
    '                    For i = 0 To Me.LoadedRawData.Matrix.i - 1
    '                        Temp5D.data(i, j, ctr2, 1) += Me.LoadedRawData.Matrix.data(i, j, k, l)
    '                    Next
    '                Next
    '                ctr2 += 1

    '            End If
    '            counter += 1
    '        Next
    '    Next

    '    'Compute perfusion values
    '    Dim PerfusionValues As New Matrix3DSingle(Me.LoadedRawData.Matrix.i, Me.LoadedRawData.Matrix.j, Me.LoadedRawData.Matrix.k)

    '    For k = 0 To PerfusionValues.z - 1
    '        For j = 0 To PerfusionValues.y - 1
    '            For i = 0 To PerfusionValues.x - 1
    '                If (Temp5D.data(i, j, k, 0) <> 0) Then
    '                    PerfusionValues.data(i, j, k) = (Math.Abs(Temp5D.data(i, j, k, 1) - Temp5D.data(i, j, k, 0))) / (Temp5D.data(i, j, k, 0) * (Me.LoadedRawData.NumberOfVolumes))
    '                Else
    '                    PerfusionValues.data(i, j, k) = 0
    '                End If
    '            Next
    '        Next
    '    Next

    '    Temp5D.Dispose()

    '    Return PerfusionValues
    'End Function

    Public Function PartitionData() As Matrix3DSingle

        'Declare Counter for RAW data
        'Declare Tagged & UnTagged Volumes
        Dim Counter As Integer = 0

        Dim TaggedCounter As Integer = 0
        Dim UnTaggedCounter As Integer = 0

        Dim TaggedVolume, UntaggedVolume As New Matrix3DSingle(Me.LoadedRawData.Matrix.i, Me.LoadedRawData.Matrix.j, CInt(Me.LoadedRawData.Matrix.k / 2))

        For l As Integer = 0 To LoadedRawData.NumberOfVolumes - 1
            'In Each Volume
            UnTaggedCounter = 0
            TaggedCounter = 0
            For k As Integer = 0 To Me.LoadedRawData.Matrix.k - 1
                'In Each Slice
                'Check if Slice is Tagged or Untagged
                If IndexArray(Counter) = 1 Then
                    For j As Integer = 0 To Me.LoadedRawData.Matrix.j - 1
                        For i As Integer = 0 To Me.LoadedRawData.Matrix.i - 1
                            UntaggedVolume.data(i, j, UnTaggedCounter) += Me.LoadedRawData.Matrix.data(i, j, k, l)
                        Next
                    Next
                    UnTaggedCounter += 1
                ElseIf IndexArray(Counter) = 4 Then
                    For j As Integer = 0 To Me.LoadedRawData.Matrix.j - 1
                        For i As Integer = 0 To Me.LoadedRawData.Matrix.i - 1
                            TaggedVolume.data(i, j, TaggedCounter) += Me.LoadedRawData.Matrix.data(i, j, k, l)
                        Next
                    Next
                    TaggedCounter += 1
                End If
                Counter = Counter + 1
            Next
        Next

        Dim perfusionvolume As New Matrix3DSingle(Me.LoadedRawData.Matrix.i, Me.LoadedRawData.Matrix.j, CInt(Me.LoadedRawData.Matrix.k / 2))

        For k As Integer = 0 To perfusionvolume.z - 1
            For j As Integer = 0 To perfusionvolume.y - 1
                For i As Integer = 0 To perfusionvolume.x - 1
                    If UntaggedVolume.data(i, j, k) <> 0 Then
                        perfusionvolume.data(i, j, k) = Math.Abs((TaggedVolume.data(i, j, k) - UntaggedVolume.data(i, j, k)) / (Me.LoadedRawData.NumberOfVolumes * UntaggedVolume.data(i, j, k)))
                    Else
                        perfusionvolume.data(i, j, k) = 0
                    End If
                Next
            Next
        Next
        Return perfusionvolume
    End Function
    Public Function RemoveNoise() As Short(,,)
        Dim reversible As Boolean
        If MsgBox("Do you want this filter to be reversible?", MsgBoxStyle.YesNo) = MsgBoxResult.Yes Then
            reversible = True
        Else
            reversible = False
        End If

        Dim threshold As Integer = CInt(CInt(InputBox.Show("Enter Minimum threshold (%)", "Minimum", "1").Text) * (LoadedRawData.Matrix.Maximum - LoadedRawData.Matrix.Minimum) / 100)
        Dim threshold2 As Integer = CInt(CInt(InputBox.Show("Enter Maximum threshold (%)", "Maximum", "50").Text) * (LoadedRawData.Matrix.Maximum - LoadedRawData.Matrix.Minimum) / 100)

        Dim temp As New Matrix4DInt16(LoadedRawData.Matrix.data)

        For i As Integer = 0 To temp.i - 1
            For j As Integer = 0 To temp.j - 1
                For k As Integer = 0 To temp.k - 1
                    For l As Integer = 0 To temp.l - 1
                        If LoadedRawData.Matrix.data(i, j, k, l) <= threshold Then
                            If Not reversible Then LoadedRawData.Matrix.data(i, j, k, l) = CShort(threshold)
                            If reversible Then temp.data(i, j, k, l) = CShort(threshold)
                        Else
                            If reversible Then temp.data(i, j, k, l) = LoadedRawData.Matrix.data(i, j, k, l)
                        End If
                        If LoadedRawData.Matrix.data(i, j, k, l) >= threshold2 Then
                            If Not reversible Then LoadedRawData.Matrix.data(i, j, k, l) = CShort(threshold2)
                            If reversible Then temp.data(i, j, k, l) = CShort(threshold2)
                        Else
                            If reversible Then temp.data(i, j, k, l) = LoadedRawData.Matrix.data(i, j, k, l)
                        End If
                    Next
                Next
            Next
        Next
        If Not reversible Then
            Return Me.GetData.Matrix.ExtractVolume(0)
        Else
            Return temp.ExtractVolume(0)
        End If
    End Function
End Class

Public Class ThresholdMaskApplicator
    Dim Threshold As Single
    Dim DataVal As Matrix3DSingle
    Dim DataInt As Matrix4DInt16
    Dim DataInt3 As Matrix3DInt16
    Dim i, j, k As Integer

    Sub New(ByVal Threshold As Single, ByVal MatrixToApplyMaskon As Matrix3DSingle, ByVal MatrixtoComputeMaskFrom As Matrix4DInt16)
        Me.Threshold = Threshold
        Me.DataInt = MatrixtoComputeMaskFrom
        Me.DataVal = New Matrix3DSingle(MatrixToApplyMaskon.data, MatrixToApplyMaskon.x, MatrixToApplyMaskon.y, MatrixToApplyMaskon.z)
    End Sub
    Sub New(ByVal MatrixToApplyMaskon As Matrix3DSingle, ByVal MatrixtoComputeMaskFrom As Matrix4DInt16)
        Me.DataInt = MatrixtoComputeMaskFrom
        Me.DataVal = MatrixToApplyMaskon
        ComputeThreshold()
    End Sub
    Sub New(ByVal MatrixToApplyMaskon As Matrix3DSingle, ByVal MatrixtoComputeMaskFrom As Matrix3DInt16)
        Me.DataInt3 = MatrixtoComputeMaskFrom
        Me.DataVal = MatrixToApplyMaskon
        ComputeThreshold3()
    End Sub
    Public Function ReturnThresholdded() As Matrix3DSingle
        For i = 0 To DataVal.x - 1
            For j = 0 To DataVal.y - 1
                For k = 0 To DataVal.z - 1
                    If DataInt.data(i, j, 2 * k, 0) < Threshold Then
                        DataVal.data(i, j, k) = 0
                    End If
                Next
            Next
        Next
        Return DataVal
    End Function
    Public Function ReturnThresholdded3() As Matrix3DSingle
        For i = 0 To DataVal.x - 1
            For j = 0 To DataVal.y - 1
                For k = 0 To DataVal.z - 1 Step 1
                    If DataInt3.data(i, j, 2 * k) < Threshold Then
                        DataVal.data(i, j, k) = 0
                    End If
                Next
            Next
        Next
        Return DataVal
    End Function

    Public Sub ComputeThreshold()
        MsgBox("Maximum Intensity Value: " & DataInt.Maximum & ". Minimum Intensity Value: " & DataInt.Minimum & ".")
        If MsgBox("Threshold using Max-Min Ratio :", MsgBoxStyle.YesNoCancel) = MsgBoxResult.Yes Then
            Dim Temp As Single = CSng(InputBox.Show("Enter Percentage: ").Text)
            Threshold = (DataInt.Maximum - DataInt.Minimum) * Temp / 100
        ElseIf CBool(MsgBoxResult.No) Then
            Threshold = CSng(InputBox.Show("Enter Threshold: ").Text)
        End If
    End Sub
    Public Sub ComputeThreshold3()
        MsgBox("Maximum Intensity Value: " & DataInt3.Maximum & ". Minimum Intensity Value: " & DataInt3.Minimum & ".")
        If MsgBox("Threshold using Max-Min Ratio :", MsgBoxStyle.YesNoCancel) = MsgBoxResult.Yes Then
            Dim Temp As Single = CSng(InputBox.Show("Enter Percentage: ").Text)
            Threshold = (DataInt3.Maximum - DataInt3.Minimum) * Temp / 100
        Else
            Threshold = CSng(InputBox.Show("Enter Threshold: ").Text)
        End If
    End Sub
End Class

Public Class HeaderOnlyExtractor
    Public ExtendedPARHeader As ExtendedPARfileInfo
    Public Sub New(ByVal FileName As String)
        Dim dti_slice_info(0) As RECimagesInfo
        Dim pixel_size As Integer
        Dim rs, ri As Single 'needed to produce display 

        ExtendedPARHeader = OpenOptimizedPAR(FileName, dti_slice_info, pixel_size, rs, ri)
    End Sub
    Function OpenOptimizedPAR(ByVal fname As String, ByRef recc_info() As RECimagesInfo, ByRef pixel_size As Integer, ByRef rs As Single, ByRef ri As Single) As ExtendedPARfileInfo
        Dim fs As FileStream 'Declare new stream

        fs = New IO.FileStream(fname.ToString, FileMode.Open)
        Dim r As StreamReader = New StreamReader(fs) 'create reader
        'read data into structure
        Dim par_struct As New ExtendedPARfileInfo
        Dim temp As String
        Dim temp2() As Single
        r.ReadLine()                            '# === DATA DESCRIPTION FILE ======================================================
        r.ReadLine()                            '#
        r.ReadLine()                            '# CAUTION - Investigational device.
        r.ReadLine()                            '# Limited by Federal Law to investigational use.
        r.ReadLine()                            '#
        temp = r.ReadLine()                     '# Dataset name: e:\intera\TH_E004_3_1
        par_struct.dataset_name = Mid(temp, 17, Len(temp) - 17 + 1)
        r.ReadLine()                            '#
        r.ReadLine()                            '# CLINICAL TRYOUT             Research image export tool     V3
        r.ReadLine()                            '#
        r.ReadLine()                            '# === GENERAL INFORMATION ========================================================
        r.ReadLine()                            '#
        temp = r.ReadLine                       '.    Patient(Name)                      :   TH_E004_11292004()
        par_struct.patient_name = Mid(temp, 45, Len(temp) - 45 + 1)
        temp = r.ReadLine                       '.    Examination(Name)                  :   Eureka()
        par_struct.examination_name = Mid(temp, 45, Len(temp) - 45 + 1)
        temp = r.ReadLine                       '.    Protocol name                      :   T1W/3D/TFE
        par_struct.protocol_name = Mid(temp, 45, Len(temp) - 45 + 1)
        temp = r.ReadLine                       '.    Examination date/time              :   2004.11.29 / 15:27:35
        par_struct.exam_time_date = Mid(temp, 45, Len(temp) - 45 + 1)
        temp = r.ReadLine                       '.    Acquisition nr                     :   3
        par_struct.acquisition_nr = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Reconstruction nr                  :   1
        par_struct.reconstruction_nr = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Scan Duration [sec]                :   113
        par_struct.scan_duration_secs = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Max. number of cardiac phases      :   1
        par_struct.max_no_of_cardic_phases = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Max. number of echoes              :   1
        par_struct.max_no_of_echoes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Max. number of slices/locations    :   60
        par_struct.max_no_of_slices = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Max. number of dynamics            :   1
        par_struct.max_no_of_dynamics = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Max. number of mixes               :   1
        par_struct.max_no_of_mixes = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Image pixel size [8 or 16 bits]    :   16
        par_struct.image_pixel_size = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        pixel_size = CInt(par_struct.image_pixel_size)
        temp = r.ReadLine                       '.    Technique                          :   T1TFE()
        par_struct.technique = Mid(temp, 45, Len(temp) - 45 + 1)
        temp = r.ReadLine                       '.    Scan mode                          :   3D
        par_struct.scanmode = Mid(temp, 45, Len(temp) - 45 + 1)
        temp = r.ReadLine                       '.    Scan resolution  (x, y)            :   240  150
        temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 2)
        par_struct.scan_resolution_x_y(0) = CInt(temp2(0))
        par_struct.scan_resolution_x_y(1) = CInt(temp2(1))
        temp = r.ReadLine                       '.    Scan percentage                    :   100
        par_struct.scan_percentage = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Recon resolution (x, y)            :   256  256
        temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 2)
        par_struct.recon_resolution_x_y(0) = CInt(temp2(0))
        par_struct.recon_resolution_x_y(1) = CInt(temp2(1))
        temp = r.ReadLine                       '.    Number of averages                 :   1
        par_struct.no_of_averages = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Repetition time [msec]             :   9.31  
        par_struct.repetition_time = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    FOV (ap,fh,rl) [mm]                :   230.00  120.00  230.00
        temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
        par_struct.fov = temp2
        temp = r.ReadLine                       '.    Slice thickness [mm]               :   2.00
        par_struct.slice_thickness = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Slice gap [mm]                     :   0.00
        par_struct.slice_gap = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Water Fat shift [pixels]           :   3.05
        par_struct.water_fat_shift = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Angulation midslice(ap,fh,rl)[degr]:   0.00  0.00  0.00
        temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
        par_struct.angulation_midslice = temp2
        temp = r.ReadLine                       '.    Off Centre midslice(ap,fh,rl) [mm] :   -0.16  45.35  0.00
        temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
        par_struct.off_center_midslice = temp2
        temp = r.ReadLine                       '.    Flow compensation <0=no 1=yes> ?   :   0
        par_struct.flow_compensation = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Presaturation     <0=no 1=yes> ?   :   1
        par_struct.pre_saturation = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Cardiac frequency                  :   0
        par_struct.cardiac_frequency = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Min. RR interval                   :   0
        par_struct.min_rr_interval = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Max. RR interval                   :   0
        par_struct.max_rr_interval = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Phase encoding velocity [cm/sec]   :   0.00  0.00  0.00
        temp2 = return_nos_from_string(Mid(temp, 45, Len(temp) - 45 + 1), 3)
        par_struct.phase_encoding_velocity = temp2
        temp = r.ReadLine                       '.    MTC               <0=no 1=yes> ?   :   0
        par_struct.mtc = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    SPIR              <0=no 1=yes> ?   :   0
        par_struct.spir = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    EPI factor        <0,1=no EPI>     :   1
        par_struct.epi_factor = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    TURBO factor      <0=no turbo>     :   150
        par_struct.turbo_factor = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Dynamic scan      <0=no 1=yes> ?   :   0
        par_struct.dynamic_scan = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Diffusion         <0=no 1=yes> ?   :   0
        par_struct.diffusion = CInt(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Diffusion echo time [msec]         :   0.00
        par_struct.diffusion_echo_time = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
        temp = r.ReadLine                       '.    Inversion delay [msec]             :   0.00
        par_struct.inversion_delay = CSng(Mid(temp, 45, Len(temp) - 45 + 1))
         Return par_struct
    End Function

End Class

Public Class ExtendedPARfileInfo
    Public dataset_name As String
    Public patient_name As String
    Public examination_name As String
    Public protocol_name As String
    Public exam_time_date As String
    Public acquisition_nr As Integer
    Public reconstruction_nr As Integer
    Public scan_duration_secs As Integer
    Public max_no_of_cardic_phases As Integer
    Public max_no_of_echoes As Integer
    Public max_no_of_slices As Integer
    Public max_no_of_dynamics As Integer
    Public max_no_of_mixes As Integer
    Public image_pixel_size As Integer
    Public technique As String
    Public scanmode As String
    Public scan_resolution_x_y(2) As Integer
    Public scan_percentage As Integer
    Public recon_resolution_x_y(2) As Integer
    Public no_of_averages As Integer
    Public repetition_time As Single
    Public fov(3) As Single
    Public slice_thickness As Single
    Public slice_gap As Single
    Public water_fat_shift As Single
    Public angulation_midslice(3) As Single
    Public off_center_midslice(3) As Single
    Public flow_compensation As Integer
    Public pre_saturation As Integer
    Public cardiac_frequency As Integer
    Public min_rr_interval As Integer
    Public max_rr_interval As Integer
    Public phase_encoding_velocity(3) As Single
    Public mtc As Integer
    Public spir As Integer
    Public epi_factor As Integer
    Public turbo_factor As Integer
    Public dynamic_scan As Integer
    Public diffusion As Integer
    Public diffusion_echo_time As Single
    Public inversion_delay As Single
    Public FOVCoordinate, Dimensions As coordinate
End Class

'Public Class BasicPixelSizeMatrixInterpolation
'    Public Function ResizePixels(ByVal OriginalMatrix As Matrix3DInt16, ByVal FinalX As Integer, ByVal FinalY As Integer, ByVal FinalZ As Integer) As Matrix3DInt16

'        Dim TempResult As New Matrix3DInt16(FinalX, FinalY, FinalZ)

'        'suppose final is 10 10 10
'        'original is 5 5 5
'        'then ratio is 2 2 2
'        'so when looking for 10 10 10 value, you are lookign for 5 5 5 value
'        'ie value for i/ratiox j / ratioy k / raioz value
'        Dim RatioX As Double = FinalX / OriginalMatrix.x
'        Dim RatioY As Double = FinalY / OriginalMatrix.y
'        Dim RatioZ As Double = FinalZ / OriginalMatrix.z

'        Dim MaxOx, MaxOy, MaxOz As Integer
'        MaxOx = OriginalMatrix.x - 1
'        MaxOy = OriginalMatrix.y - 1
'        MaxOz = OriginalMatrix.z - 1

'        Dim temp1, temp2, temp3 As Integer
'        For i As Integer = 0 To FinalX - 1
'            For j As Integer = 0 To FinalY - 1
'                For k As Integer = 0 To FinalZ - 1

'                    'check if the ratio isnt out of the original matrix
'                    temp1 = CInt(IIf(i / RatioX < MaxOx, i / RatioX, MaxOx))
'                    temp2 = CInt(IIf(j / RatioY < MaxOy, j / RatioY, MaxOy))
'                    temp3 = CInt(IIf(k / RatioZ < MaxOz, k / RatioZ, MaxOz))

'                    TempResult(i, j, k) = OriginalMatrix(temp1, temp2, temp3)
'                Next
'            Next
'        Next

'        Return TempResult
'    End Function
'End Class

